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A study based on inviscid analysis has been made to examine the 
flow field produced from a convergent^divergent nozzle when a. strong 
curved shock occurs. It is learned that a certain constraint is im- 
posed on the flow solution of the problem which is the unique feature 
of the flow within this flow regime and provides the reason why the 
inverse method of calculation cannot be employed for these problems. 
An approximate method has been developed to calculate the flow field 
and results have been obtained for two-dimensional flows. Analysis 
and calculations have also been performed for flows with axial sym- 
metry. It is learned that under certain conditions, the vorticity 
generated at the jet boundary may become infinite and the viscous 
effect becomes important. Under other conditions, the asymptotic 
free jet height as well as the corresponding shock geometry have 
been determined. 
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Subscripts 


••• ■ NOMENCLATURE 
coefficients of <J> function 
Mach number 

curvilinear coordinates measured normal to and along streamline 

pressure 

velocity 

gas constant 

entropy 

time 

velocity components in x~ and y -direct ions 

Cartesian coordinates 

shock profile 

shock turning angle 

dimensionless normal coordinate 

stream angle 

ratio of specifi c heats, ic = 1,4 for air 

dimensionless horizontal coordinate 

density 

shock angle 

velocity profile 

stream function 

vorticity 


a 

c 
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asymptotic state conditions 
evaluated at nozzle centerline 
evaluated at nozzle exit plane 
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j evaluated at free j et boundary 

o stagnation conditions 

s evaluated at curved shock 
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1. INTRODUCTION 
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The study of compressible flow through a convergent -divergent nozzle 
provides the most basic and fundamental knowledge of the flow of a com- 
pressible fluid. Depending upon the ambient pressure ratio (ambient to 
the supply pressure), a variety of flow patterns occurs both, inside and 
outside the nozzle. While the flow inside the nozzle will be essentially 
unchanged for sufficiently low pressure ratios, the free jet flow patterns 
may be entirely different due to the existence of a family of oblique shock 
solutions. These flow phenomena may be further complicated by the fact 
that the real fluid .is viscous and the flow may separate (shock induced 
separation) from the solid wall inside the nozzle when the wall boundary 
layer cannot cope with the high ambient pressure prevailing outside the 
nozzle. One should recognize, however, that in order to understand the 
viscous effect on the overall flow pattern, it is a prerequisite that the 
flow solution based on purely inviscid analysis be established. 

In the major operating range of the nozzle pressure ratios, the in- 
viscid flow patterns of the free jet flow are relatively simple, and the 
description of them can be found in standard text, books [1,2,3].* There 
exists, however, a range of the pressure ratios where the inviscid free 
jet flow patterns are very complicated and have not been established. 

Referring to Fig. 1 where a series of two-dimensional nozzle free jet 
flow is depicted, an oblique shock is generated from the edge of the noz- 
zle when the ambient pressure p is kept slightly above the exit pressure 
p . This shock can be regularly reflected from the centerline (Fig, la) 
if the strength of the incident shock is not strong enough to prevent the 
occurrence of such a reflection. For higher pressure ratios where such a 


^Numbers in brackets refer to entries in REFERENCES. 
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regular reflection is not possible, Mach reflection may occur within 
the flow field (Fig. lb) . Recent studies by Chang [4] and Chow and 
Chang [5,6] on Mach reflection of shock and Mach disc showed that for 
small nozzle Mach numbers (M^ < 1.484 for k = 1.4), Mach reflection 
cannot occur since the triple point condition based on inviscid analy- 
sis cannot be satisfied. For larger nozzle Mach numbers, Mach reflec- 
tion of shock is possible and the Mach stem height will be increased 
for higher pressure ratios. Even in this situation, there exists a 
range of nozzle Mach numbers where the Mach stem height can never reach 
the height of the nozzle for sufficiently high pressure ratios and a 
strong curved shock necessarily occurs in the flow field (Fig. 1c) . 

Perhaps the above description of the flow events may be best illus- 
trated in Figs. 2, 3, and 4. Figure 2 shows the possible conditions for 
the occurrence of a triple point (Mach reflection pattern) and various 
limits for shock reflection. It is also recognized that from inviscid 
considerations no triple point can occur for M < 1.484 (k = 1.4). 

Figure 3 illustrates the Mach stem height calculated by either detailed 
calculations or by approximate method corresponding to different operat- 
ing conditions. It may be seen that for moderate supersonic nozzle Mach 
numbers, the maximum Mach stem height is still smaller than the height 
of the nozzle. Figure 4 illustrates the corresponding pressure ratio when 
Mach reflection occurs. It is now obvious that for pressure ratios higher 
than those required for Mach reflection (or regular reflection for M < 
1.65) but lower than the corresponding normal shock values, the regime 
of strong curved shock prevails.* Under this condition, a shock wave 


*This regime may overlap that of Mach reflection. 
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(which may even belong initially to the weak shock solution) is generated 
at the corner of the nozzle (Fig. 1c) and proceeds downstream with increas- 
ing wave angle and shock strength until the normal shock condition is at- 
tained on the centerline of the nozzle. This shock is never reflected, 
and the downstream flow field is bounded by a constant pressure boundary. 
The majority’ of the flow field is subsonic. Since the curved shock gener- 
ates vorticity, the flow is thus rotational. Eventually all streamtubes 
are straight and parallel to the nozzle centerline so that the pressure 
becomes uniform and is equal to the ambient, pressure. However, the ve- 
locity, density, etc., are non-uniform which is, of course, the result of 
the existence of the vorticity. This state-- termed the asymptotic state-- 
can exist only at the mathematically ’’far downstream" positions. 

It should be pointed out that in previous studies [4,5,6] on the oc- 
currence of triple point, inviscid considerations were found to be ade- 
quate from experimental observations made with a shock tube [7] for rela- 
tive strong incident shocks.. However, for weak incident shocks (small 
Mach number and small shock strength), triple point has been observed [7] 
under conditions even when the inviscid analysis does not yield a solution. 
It is conceivable that under these conditions, the thickness of the shock 
wave may not be thin enough that its subsequent interaction with the flow 
fields may become an important point of consideration for the solution at 
the triple point. In view of these observations [7], the range of occur- 
rence of Mach reflection presented by Chow and Chang may not be precise; 
nevertheless, it provides useful guidance in describing the flow regimes 
of free jet flow under different operating conditions. 

While the existence of such a strong curved shock flow regime has 
been mentioned only briefly by Ferri [2], no work has been performed in 
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the study of associated free jet pattern? within this flow regime. It is 
true that all nozzles do not normally operate • under such flow conditions 
and even they do, the viscous effects may significantly modify the flow 
pattern. The motivation of studying such flow problems comes not only 
from filling the gap of information which, so far, is lacking. The unique 
character of such a rotational subsonic or mixed flow field downstream of 
a curved strong shock and extending to far downstream positions offers con- 
siderable interest and challenge. It should be reiterated that a flow so- 
lution, on the basis of purely inviscid flow analyses, should be estab- 
lished before any viscous effects can be incorporated into the consider- 
ations. This is indeed the purpose of the present investigation. 

Since the problem under consideration concerns the determination of 
the profile of a strong curved shock and the details of the free jet flow 
behind it, it is similar to the blunt body problem, i.e., a blunt body 
immersed in a supersonic flow with a detached shock wave. The method of 
analysis for the blunt body problem will be reviewed. 

The blunt body problem has been a matter of concern since early 1950, 
and many methods for solving problems of this type have been developed. 
These methods of approach may be classified as solutions to direct or in- 
verse problems. For the direct problem, the shape of the body is given 
while the shock profile and details of the flow field between the shock 
and body are to be examined. In the inverse problem, the shock shape is 
specified and the resulting body configuration and accompanying velocity 
and pressure fields on, or ahead of, the body are to be calculated. These 
different methods of analysis are discussed in the following section. 

For the direct problem, the analyses usually can be divided into 
(1) the method of integral relations, (2) the unsteady flow approach, and 
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(3) the relaxation method. "'The method of integral relations was originally 
developed by Dorodnitsyn [8] as a general method of numerical solution for 
nonlinear hydrodynamic problems. It reduces the problem of integrating a 
system of nonlinear partial differential equations to that of solving an 
approximate system of ordinary differential equations. Phis metnod was 
originally applied to the blunt body problem by Belotserkovskii [9] to. 
calculate the supersonic flow of a perfect gas past a circular cylinder. 
Following this method of approach, many calculations for more general body 
shapes in symmetric flows have also appeared [10] . 

The use of the method of integral relations was also proposed by 
Mel’mikov [11] for calculating the flow behind the Mach stem of an over- 
expanded supersonic nozzle where Mach reflection occurs. For the case of 
a single strip, a system of ordinary differential equations with corre- 
sponding boundary conditions describing the problem was derived. However, 
no numerical results were obtained. 

In the unsteady flow method,, the steady flow solution is obtained as 
the limiting state is reached asymptotically in time from an unsteady flow 
with constant free stream and body conditions and suitably chosen initial 
conditions. The unsteady flow equations are solved by finite difference 
techniques. Von Neumann and Richtmeyer [12] first suggested this idea to 
calculate one-dimensional flow associated with shocks. In their calcu- 
lation, large computer storage is required and, as a result of the shock 
not being considered as a sharp discontinuity, the computational time is 
too lengthy. Later, Von Neumann and Richtmeyer r s idea was developed and 
extended by Lax and Wendroff [13] and others. By properly setting the 
equations into finite difference form, soluti.ons may be obtained in which 
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shock waves are approximated by narrow regions where the physical parame- 
ters undergo drastic changes. Moretti and Abbett [14] further improved 
this method by introducing special three-dimensional method .of character- 
istics to describe the conditions on the shock and body boundaries. The 
points between shock and body were calculated by employing the technique 
developed by Lax and IVendroff. This method has also been extended to 
calculate the three-dimensional steady flows by Moretti and Bleich [15]. 

In the relaxation method, a stream function is usually introduced 
as a dependent variable. The differential equation for the stream func- 
tion is considered as the basic equation. The plane for numerical cal- 
culations is covered by a network of squares. The basic differential 
equation is written in finite difference form and a residue which must 
be zero everywhere when the final solution is reached is defined at each 
point. The assumed values of the dependent variables are then substituted 
into these equations. In general, the residue which is a measure of the 
failure of the assumed solution to satisfy the equation will not be zero. 
The problem is to reduce by suitably adjusting the assumed values of the 
dependent variables the residuals to zero at every point. Gravalos, 
Edelfelt, and Emmons [16] and Hayes and Probstein [17] used this method 
to calculate the blunt body problem. 

In the inverse problem, the shape of the detached shock is assumed 
and the numerical calculations are subsequently carried out until the 
body whose geometry is obtained by locating the points corresponding to 
constant stream functional values (e.g. , \J) = 0) is reached. One may ob- 
serve that the resulting body shape corresponds to the selected detached 
shock geometry. However, fundamental questions arise with respect to the 
uniqueness and the existence of a solution with respect to the stability 
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convergence of calculation, procedures. Hayes and Probstein [17} stated: 
"the most important single feature of this inverse problem is 
the insensitiveness of the shock shape to local changes in 
the body shape, and this feature leads to essential diffi- 
culties in the inverse problem. A minute local change in the 
shock shape will generally cause a large change in the body 
shape and may even preclude the existence of a solution." 

From a mathematical point of view, they also stated: 

"the determination of the subsonic part of the flow field in 
the problem is governed by an elliptic partial differential 
equation" 

and 

"we know as far as elliptic differential equations are con- 
cerned, the initial value problem is improperly posed and 
leads to an unstable solution when treated by finite difference." 
However, the inverse problem was solved successfully by several authors 
[18,19,20] who found that suitable numerical methods would yield suf-^ 
ficiently accurate results even though such methods have regions of 
numerical instability. 

It is now appropriate to point out that for blunt body problems, 
the elliptic region is finite. The flow will pass through a sonic line 
and become supersonic somewhere along the body. This may be the reason 
that, even when numerical instability and insensitivity of shock to body 
shapes exit, indirect problems can be calculated successfully. For the 
present problem, the elliptic region extends to infinity and is thus un- 
bounded. It will be seen that application of such an inverse method, 
i.e., guessing a shock shape and "marching downstream," would be completely 
unworkable; in fact, it is doomed to fail. 
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Since the present flow regime of curved strong shock is identified 
only after the phenomenon of Mach reflection is understood, it is also 
necessary to describe briefly the work by Chow and Chang [6] on Mach 
reflection. It was recognized that when Mach reflection occurs, the 
central core flow will pass through a strong shock (Mach shock) while 
the fluid above it will go through relatively weak incident and re- 
flected shocks. Since it is necessary that these two portions of the 
fluid assume the same pressure and same streamline flow angle along 
the common boundary, it is obvious that the mutual interaction between 
the streams plays the controlling role to the solution of the problem. 

It was shown that for supersonic flow state behind the reflected shock, 
the core flow will eventually pass through a throat-like condition simi- 
lar to the choked secondary flow within the supersonic ejector system [21] 
The detailed calculations of the flow field, including the determination 
of the Mach stem height, may be performed with the method of characteris- 
tics for the upper stream and with the simplified one-dimensional analysis 
for the core flow after the method of integral relations has been applied 
for both streams. However, when the state behind the reflected shock is - 
subsonic, such detailed calculations are not possible and the Mach stem 
is estimated from an approximate method stressing the interacting force 
between two streams. These considerations have also been extended into 
axisymmetric flow configurations. Some of these results are presented in 
Figs. 2, 3, and 4; other results were reported by Chang [4] and by Chow 
and Chang [6]. 
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2 . THEORETICAL CONSIDERATIONS 

2.1 DESCRIPTION OF THE PROBLEM 

On the basis of the discussion presented in the preceding chapter, 
it is recognized that the present problem occurs for a supersonic flow 
issuing from a nozzle when the ambient pressure is high enough that an 
oblique shock is generated at the corner of the nozzle. The wave angle 
of the shock increases as it moves downstream until a normal shock is 
attained on the centerline of the nozzle. The flow behind the shock is 
bounded on one side by a constant pressure free jet boundary and is ro- 
tational. The viscous effects occur ing along the jet boundary and within 
the fluid are disregarded. Eventually, all streamlines become parallel 
to the axis of symmetry, and the pressure is equal to the ambient value. 
The velocity, however, is non-uniform at this asymptotic state as a re- 
sult of the presence of vortici ty. A schematic diagram depicting such a 
flow field is shown in Fig. 5. Since the flow is symmetric with respect 
to the x-axis, only one-half of this region will be considered. 

2.2 BASIC RELATIONS GOVERNING THE FLOW 

One now proceeds to formulate the problem in a more detailed and 
precise manner. The governing relations are listed sequentially in 
the following sections. 

2.2.1 Shock Relations 

For the occurrence of an oblique shock within a supersonic 
flow field, the governing algebraic equations are: 
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2.2.2 Free Jet Flow Region 

For a two-dimensional invisciu flow of a compressible fluid., 
the continuity and momentum principles can be given as: 

(Hereafter, unless mentioned otherwise, all fluid properties will be 
treated as normalized by the corresponding values at the nozzle exit 
condition, and all length quantities are normalized by the half-height 
of the nozzle.) 
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The energy equation can be given simply as 
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where S denotes entropy. It is understood that due to the curved shock 
prevailing upstream, the flow field is not homentropic. Nevertheless, 
the flow is isentropic along individual streamlines. This would also 
imply that the flow is isoenergetic throughout the field. 

The flow region under consideration is bounded by (1) a constant 
pressure jet boundary whose profile, y^ (x) is yet unknown, (2) a line of 
symmetry which is the x-axis, (3) the upstream curved shock whose shape, 
x (y) , is also unknown, and (4) the asymptotic state prevailing at far 
downstream position (i.e., x -*• “) . The boundary conditions imposed on 
the solution of Eqs. (2.6), (2.7), and (2.8) are: 


p p a along y. (x) 

u = u 'h 

s 

V = V 

along x s (y) 

P = P s I 



( 2 . 10 ) 


( 2 . 11 ) 


and 

v = 0, ~ = 0 at y = 0 (2.12) 

ay - 

In addition, any solution corresponding to a particular ambient pres- 
sure ratio must be subjected to a certain constraining condition when the 
asymptotic state is reached. This is considered in detail in the next 
section. 


2.2.3 Constraining Condition at the Asymptotic State 

One now gives attention to the control volume shown in Fig. 6a. 


The continuity principle would yield 
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The x- momentum principle also gives 
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It should be noted that for the present problem, the right-hand side of 
Eq. (2.15) is a constant value which is greater than zero but less than 
unity . 

After introducing p . , u . , and y . as reference quantities for the 

aj aj aj 

integrands in Eqs . (2.13) and (2.15) and defining = u a / u a j an d 

q = y/y these equations become 
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Upon combining Eqs. (2.16) and (2.17), one obtains 
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and C aj is a constant for a particular nozzle Mach number operating under 

a particular pressure ratio. It should be noted that the two integrals 

have positive definite integrands which are never greater than unity as 

(j) varies from <f> (<j> = u /u . < 1) to unity: thus, 
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Equation (2.19) provides the constraining condition of the flow when the 
asymptotic state is reached which in turn imposes the additional condition 
.which the solution of the problem must satisfy. 

It is now clear that this constraint is the reason that the indirect 
approach to the blunt body problem cannot be adopted for the present 
calculations. The hope of estimating a curved shock profile and marching 


downstream and correcting the upstream shock shape from some unmatched 
properties prevailing downstream can never be realized-. If a certain shock 
shape is assumed, it immediately determines the corresponding asymptotic 
state for each of the stream tubes. Usually this profile will not satisfy 
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the constraining condition 'given by Eq. (2.19). Marching downstream 
with such a shock will necessarily result in illogical mathematic 
operations (e.g., seeking the square root of a negative number) if these 
calculations were performed on a digital computer and the computation 
was interrupted before the asymptotic state was reached. Since an ac- 
curate estimation of the shock shape* is impossible, such interruption 
can never be avoided. This explains why the indirect method cannot be 
applied to the present problem. 

A method has been devised and is suggested to solve the present 
flow problem. It is recognized that with a given nozzle operating at 
a particular pressure ratio within this flow regime, the asymptotic* con- 
ditions of the jet boundary and center streamlines are known. A velocity 
profile at the asymptotic state which, in addition to fulfilling all re- 
quired boundary conditions, contains one additional parameter which is to 
be determined by the constraining relations given by Eq. (2.19). Once 

the asymptotic profile and thus y . are determined, the associated up- 

aj 

stream curved shock geometry can be evaluated and the flow field between 
the shock and the asymptotic state can be determined through numerical 
calculations. 

i 

2.3 ASYMPTOTIC STATE 

Since the solution of the problem relies on the establishment of the 
asymptotic flow profile, a detailed examination of the conditions imposed 
on such an asymptotic flow state is therefore necessary. 

*It is believed that even if the correct shock shape were inserted upstream, 
the subsequent error of calculations due to round-off or truncation would 
eventually cause the termination of calculations. 
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Under the asymptotic, .flow condition; the slope of the velocity pro- 
file is precisely the vorticity of the flow. A study of the vorticity 
generated by the shock* at the corner of the nozzle would also yield 
the slope of the profile of the jet boundary streamline at the asymptotic 
state . 


2.3.1 Vorticity at the Jet Boundary Streamline 

One now examines the flow condition (downstream of the shock) 
at the corner of the nozzle where the shock is generated (Fig- 7) . In 
the streamline coordinate system, the continuity equation is given in 
dimensional form by: 


1 - M a 

6 la + £0 = 0) 

q 3s 3n 


( 2 , 21 ) 


and the equations of motion for the s and n directions are, respectively: 

n M = _ I (2.22). 


3s 


and 


p q 


2 

3s' 


P 3s' 


3p 
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(2,23) 


where 0 is the streamline angle and q is the velocity of the fluid. The 
vorticity can also be expressed by 


_ 3q 30 _ p 3S t 
W 3n ^ 3s p q R 3n' 

For a j et boundary streamline at the corner of the nozzle, Eq , 
implies 30/3nk = 0 as 3p/3s = 0. It may be seen that 


(2.24) 


( 2 . 21 ) 


„ = IS = P / dS da 30 \ = Q 

j p q R 3n p q R ' da d0 Bn * 


(2,25) 


*The infinite vorticity due to the discontinuous change in velocity at the 
jet boundary is, of course, not under consideration. 
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It is thus recognized that ^vanishing vorticity occurs along the jet 
boundary, implying that the slope of the velocity profile at the asymptotic 
state should vanish at the jet boundary. 


2.3.2 Velocity Profile at the Asymptotic State 

One now realizes that this velocity profile cf> in Eq. (2,19) 
should satisfy the conditions 


(.i) At n = 0: 


(ii) At n = 1: 


* = +c 


u 

ac 


u . 


ai 


and 


dn 


o 


(j) = 1 and - 0 


(2,26) 


(2. >27) 


In view of the simple requirements to be satisfied by such a profile, 

many functions seem to satisfy the need of the present analysis. It should 

be noted that it is not a simple matter to select such a profile. It is 

pointed out that a very stringent requirement which has not been presented 

2 2 

herein is that the curvature (or d“<|>/dq ) of the profile should be non- 
negative at ri = 0 for the whole range of the pressure ratios. As no 
physically possible flow situation can occur under the. condition when 
d <f>/dp is negative, it has been found that polynomials, sine, and co- 

sine functions , cannot meet this requirement throughout the range of pres- 


sure ratios. 

A profile has been selected which has the following form: 

4> = A + B erf [c (1 - p 2 ) ] + D erf [c (1 - p 4 ) ] (2.28) 

where A, B, and D are parameters to be determined from conditions given by 
Eqs . (2.26) and (2.27). Only even powers of n appear in the above equation 
due to the condition of symmetry. Upon evaluating these parameters, 

Eqs . (2 . 28) becomes : 
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2 (1 - <y : 2 1 - V 4 

» = 1 - - e r-r cc f - erf [c(1 - n 3] erf [c(1 ' n )] (2 - 29) 

The coefficient C is yet to be determined for each particular pressure 
ratio so that the constraining relationship, Eq. (2.19), can be satis- 
fied. A range of values of C has been found for different, nozzle Mach 
numbers under various pressure ratios and is given in Table 1. 

Once the value of C is determined from the constraining condition, 

the free jet height y . at this position can also be computed and the 

aj 

detailed shape of ■ the asymptotic profile completely established. This 
information would readily lead to the determination of the geometry 
and the strength of the upstream curved shock wave. Some of these typi- 
cal velocity profiles are presented in Fig. 9. 

2.4 CORRESPONDING UPSTREAM SHOCK PROFILE 

After the asymptotic height of the free jet, y and its detailed 

a J 

profile are established, the corresponding curved shock configuration 
can be determined through the following considerations: (Also, see 

Fig. 6b.) If one carries out the integration 

U 

J p a u a dy (2.30) 

O 

along the asymptotic profile, the continuity principle would imply that 
the fluid assuming the position y at the asymptotic section would have 
passed through the shock at height y which is equal to 

P a u a dy (2.31) 

o 

since the velocity of the fluid at y has a known value (j) , its stagnation 

a a 
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pressure p can be found -from 
oa 


( k/ tc- 1) 


oa 


Pa 


J. 


2 2 
1 - C . <J> 
aj 


(2.32) 


The ratio of the stagnation pressures across the shock at y can be re- 

s 

lated through the shock wave angle a by 

k + i 2 .2 

/ — r -“e sin 0 \ 


K + 1 ..2 . 2 

1 + — - — M sin a 
2 e 


p P P 
1 oa 1 oa r a 


2k ,,2 . 2 K - 

— M sm a - — 

k + 1 e K 


H) 


I/k- 1 p p p 

oe r a oe 


(2.33) 


where P a /p Qe is the operating pressure ratio. 

Thus, the wave angle a can be found as a function of any position y c 
along the shock. Additional integration of the relationship 

dy„ 


dx 


tan O 


(2.34) 


with the initial condition of = 0, y g = 1 would produce the shock con- 
figuration x g (y ) corresponding to the established asymptotic profile. 


2.5 FLOW FIELD BETWEEN THE SHOCK AND THE ASYMPTOTIC STATE 

One now proceeds to investigate the detailed flow field between the 
curved shock and the final asymptotic state. The basic differential 
equations governing the flow are rewritten as: 

3u . Bv \ Bp 


/ ou dv \ dp Bp 
p ( 157 + 77 ) * u 77 + v ^ = 0 


P (u 


Bx By 
3u 


By 


Bx + V By 


= . iE 

3v / . -2 Bx 


k M 


( 2 . 6 ) 

C2.7) 
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^ r 3v 3v 

P ( U 7 + V 7 


K M 


_ 32. 

2 8y 


( 2 . 8 ) 


For the convenience of numerical calculations, another pair of inde- 
pendent variables, E, and ip, which are related with x and y coordinates 
through 


i = 


X - x s OfO 
1 + X. - x s OJO 


dip . 3^ 

TZ = " pvand 7 = pu 


(2.35) 


(2.36) 


are introduced. i> is the dimensionless stream function and x , which can 
easily be interpreted as a function of ijj, is the upstream curved shock 
configuration. Upon adopting this pair of independent variables, the 
semi-infinite physical region is transformed into a rectangular finite 
domain in the £,\p plane. Lines of constant ip values correspond to stream- 
lines, shock corresponds to £ = 0, and ip ~ 1 corresponds to the free jet 
boundary. The asymptotic state assumes the location £ - 1 (see Fig. 8). 

It may easily be established that the transformational relations are: 

dx 

dijJ / 3£ dip 


3_ 

dx 


y 


n in 2 / ‘ UA s) 3 3 

= 0-?) p - PV — J " P v 


(2.37) 


and 


3_ 

9y 


= (1 - Q Pu 


dx ^ ^ 

s 3 3 

dip 3£ + PU dip 


(2.38) 


and Eqs, (2.6), (2.7), and (2.8) become 

dx c 


(1 - 5) 2 P + pv )f| - p 2 V §| - (l - O 2 p“ u -gf- 


t P 2 U + (1 . n 2 u ie. = 0 
dip U ^ 3^ 


(2.39) 
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... r-,2 3u 1 

(.1 - O pu ^ + 2 

^ K M 


(1-01 + pv 


dXo 


_s\ 3£ 1 nv 3£ 

df/35 ' K M 2 ^ 


9 pv = 0 


(2 . 40) 


„ n 2 3v 1 ,. 2 dx s3p 1 3p 

0 ‘ ^ 35 . 2 1:1 “ ^ dil) 35 .,2 3i() 0 

s K M + v K M 


(2.41) 


In addition, the isoenergetic condition is given by 

2 2 
v 

2 “ 


I _1_ £ u" + v“ 1 


K - 1 .,2 p 

M 

e 




1 1 


2 tC ~ 1 y.2. 

M 

e 


(2.42) 


The isentropic flow relationship along individual streamlines can be 


written as 

^ = f(i|i) C 2 - 43 ) 

P 

where f (ip) is evaluated according to the already established shock con- 
ditions . 

It may be observed that the asymptotic flow" condition which may be 
given as 

-n = r> v — D —2- — 0 — ^ — z: 0 

P p a J ’ 3£ u ’ 3 \Jj u 

|| = 0, || = 0, and |^ = 0 at 5 = 1 (2-44) 

satisfies Eqs . (2.39), (2.40), and (2.41) identically. 

The boundary conditions for the system of equations are : 


(1) 

At £ = 0, all flow properties 

are provided from the shock 



relations given by Eqs. (2.1) 

through (2.5). 

(2.45) 

(2) 

At £ = 1, all flow properties 

are provided from the 



asymptotic state. 


(2.46) 

(3) 

At t|) = 1, p = p , and p = p . 

a aj 


(2.47) 

(4) 

At ip = 0, v = 0, 3p/3ip - 0, Su/dilt = 0, and Bp/34; = 0. 

(2.48) 


ORIGINAL PAGE IS 
OP POOR QUALITY 
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2.5.1 Numerical Calculations 

A predictor-corrector finite difference method has been ap- 
plied to solve numerically the above partial differential equations (i.e., 
Eqs. (2.39), (2.40), and (2.41)) with the corresponding boundary conditions 
given by Eqs. (2.45), (2.46), (2.47), and (2.48). In order to carry out 
the numerical calculations, the square ABCD in the plane has been 
divided into small meshes as shown in Fig. 8c. 

In the predictor phase of the numerical method, flow properties p, 
u, and v within the region are estimated. The pressure can always be re- 
lated with these p, u, and v values through the isoenergetic condition or 
the isentropic relationship when its ^ value is known. To obtain the 
initial information of p, u, and v, their values along the £ = 0.5 line 
is established initially by obtaining: 

(1) The p^, u r> and v * values from the standard three-point 
central difference scheme from the known values of properties 
at B, - 0 and £ = 1, 

(2) The p , u,, and v. values from the basic equations, and 

w V ^ ^ 

(3) Integration of these derivatives from ^ = 1 toward ^ = 0 
numerically by applying the formula of 


Sf 

f 0J> - Aijo = f 010 - (ifO Ai/> 

where f denotes any representative property of the fluid. 
Initial values of u (or v) at i]j = 1 must be adjusted so that the con- 
dition 


(2.49) 


v (^0) = 0 


(2.50) 


is satisfied on the axis. 


'^Subscripts now indicate differentiation. 
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These processes are, repeated for constant E, lines between 
0 < 5 < 0.5 and 0.5 < £ < 1.0. 

In the corrector phase of the calculations, values of p, u, and v 

are recalculated by following the same procedures but using the standard 

five-point central difference formula [22] in £ derivatives and the three - 

point central difference in \p derivatives* starting from the E, = 0,5 line. 

In correcting the property values from the ^ = 1 line toward the axis 

(i|j = 0} , a more sophisticated local predictor-corrector technique [20] has 

been applied. These correcting calculations should be repeated until the 

differences of the fluid property between the successively determined 

values of all node points in the plane are within a desirable limit. 

In this study, this limit is taken as 10 

Once the final flow properties are determined, including those along 

the jet boundary, the profile of the jet may be established by integrating 

the expression dy./dx. = v./u. with the initial condition of y . - 1 at 
r J J 3 J 3 

x — 0 . 

2.6 DEGENERATE STATE 

It should be pointed out that for each nozzle Mach number, there 
exists a pressure ratio corresponding to the normal shock appearing at 
the exit section of the nozzle. The downstream flow field is uniform 
(i.e., y. = 1) and the constraining condition becomes the Rankine-Hugoniot 
relationship which is automatically satisfied by the normal shock. Such 
a state is termed the "degenerate state" and the corresponding pressure 
ratio for its appearance is presented in Fig. 4. 


*This must be modified for the lines adjacent to shock, the asymptotic 
state, and the iJj = 1 and t/; = 0 streamlines. 
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3 ... RESULTS AND DISCUSSION 

The method described in the previous chapter has been translated 
into a FORTRAN IV computer program with double precision. For a given 
initial Mach number at a certain operating pressure ratio, it takes 
approximately five seconds of computational time on the IBM 360/75 
digital computer to calculate all flow properties related to the prob- 
lem. In this study, cases of three initial Mach numbers (M = 1.49, 3.0, 
and 5.0) operating at different pressure ratios have been calculated. 

The results are plotted and presented in Figs. 9 through 18. 

In Fig. 9, the calculated asymptotic velocity profiles $ are plotted 
against the coordinate r|. It can be seen from those profiles that the 
change of <f> from <j) to unity occurs mainly in the upper part of x\. This 
is especially true for higher pressure ratios and also true for higher 
nozzle Mach numbers. From Table 1 it can be seen that the coefficient 
C of <j> function (Eq. (2.29)) has larger values at higher pressure ratios. 
This explains why the core flow is essentially uniform under these flow 
conditions . 

The results of the free jet flow calculated by- predictor-corrector 
scheme also shows that when the operating pressure ratio is sufficiently 
high, there is a uniform core flow around the centerline. The portion 
of this uniform core flow increases as the pressure ratio is increased.. 
When the normal shock pressure ratio is reached, the entire free jet 
flow becomes uniform. 

The calculated centerline pressure variation and the free jet 
boundary profile are plotted against the x coordinate in Figs. 10 and 11. 
As the operating pressure ratio increases toward normal shock pressure 
ratio, the following phenomena have been observed: 
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(1) The value of shock stand-off distance x , shown in Fig. 12, 

sc 

is reduced toward zero as the normal shock condition is ap- 
proached. 

4 (2) The corresponding asymptotic free jet boundary height, y 

a J 

increases toward unity. 

(3) The practical distance, x , for the free jet flow to reach the 

a 

asymptotic flow conditions is also reduced. 

The pressure distribution and its variation throughout the field 
corresponding to specific flow conditions are shown in Fig. 13. Example 
to illustrate the differences in flow properties between the initial pre- 
dicted value and final corrected values at £ = 0.5 are presented in 
Figs. 14 and 15. 

It should be mentioned that the asymptotic velocity profile repre- 
sented by equation 4> is certainly not the only type which can be adopted 
for the present problem. In a different effort, a profile of the form 
given by 


2 U - V 2 1 - <f> c 4 

+1 " 1 - tanh (c) tanh t c a - 1 )] + tanh (c) tanh t c (1 * 1 )] 

(3.1) 

has been employed, and the results are shown in Fig. 16. It can be seen 
that the difference between the 4> and <p profiles under the same con- 
dition is indeed small, and this small difference has also been observed 
under various pressure ratios. However, if one adopts a polynomial pro- 
file of the form 

<}> 2 = A + B q 2 + C n 4 + D n 6 , (3.2) 

it is found that the behavior of this profile is quite different from that 
of $ or 4^ at the same pressure ratios. In fact, when a certain level of 
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pressure is reached, the profile has an absolute minimum velocity ocur- 
ring not at the centerline which does not correspond to any physically 
realistic flow condition (see Fig. 17). This is why the profile of a 
simpler form given by equation (j)^ has not been employed for the present 
s tudy . 

The solution of the problem, nevertheless, relies on the selection 
of the type of asymptotic profile. Although profiles of a similar 
nature may be adopted for this study, and the final results are not much 
different from each other, these solutions can be considered as only ap- 
proximate solutions. Meanwhile, the selection of the profile in its 
present form, which is monotonical ly increasing for increasing y values, 
has dismissed the possibility of profiles with negative slope (positive 
vorticity) away from the centerline, whose occurrence may seem to be im- 
probable but cannot be ruled out by any basic principles. In this sense, 

the selection of the asymptotic profile in its present form is somewhat 

<,* 

restrictive . 

Perhaps it is due to the restrictive character of the profile that 
some numerical instability of the results has been observed. In numeri- . 
cal calculations by finite difference schemes, it is generally recognized 
that for better definition or resolution of the flow field, finer grid 
should be employed. In the present study, a small step of q (Ary = 0.01) 
has been employed for calculations since it is needed especially in the 
region close to the centerline where normal shock occurs. However, all 
results reported here were obtained with A£ = 0.125. Finer grid size in 
the £ direction has led to oscillatory flow results, Figure 18 shows one 
set of these results for the pressure distribution along the centerline 
when a smaller A£ (A£ - 0.0625) is employed. It. is thus obvious that with 
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the saute numerical limit of 10 ^ as the margin for convergence, a grid 
finer than = 0.0625 may not even lead to convergent solutions. 

It should also be mentioned that the question of how this curved 
shock flow regime merges with the regime of Mach reflection cannot be 
answered in this investigation. Additional study or even experimental 
explorations are needed before the problem can be resolved successfully. 

Similar considerations can be applied to axi symmetric flow con- 
figurations within this curved shock flow regime. A detailed study and 
results are reported in the appendix. It should be noted that with the 
axisymmetric geometry, the vorticity along the jet boundary is no longer 
zero. In fact, on some occasions, the vorticity may approach infinity. 
Under this condition, the viscous effect, of course, becomes important. 

Prompted by the occurrence of such unusual occasions, perhaps it 
should be suggested that the vorticity consideration on the jet boundary 
is not important and should not be taken into consideration in the es- 
tablishment of the asymptotic profile since viscous mixing does occur 
along the free jet boundary in actual flows. Whether this change in 
approach would result in simpler analysis or simpler asymptotic profiles 
is yet to be examined. 



4 CONCLUSIONS 


From this study, it is recognized that a strong curved shock flow 
regime exists in the operation of the convergent -divergent nozzle. Al- 
though the precise manner of change-over between the Mach reflection 
and the present curved shock flow regimes has not been clearly defined, 
the pattern of the flow field associated with the occurrence of a 
strong curved shock is, nevertheless, established, and the gap of in- 
formation on nozzles operating under these conditions is satisfac*- 
torily filled. The inviscid flow analysis of nozzle flow problems thus 
yields consistent flow patterns throughout the entire range of the pres 
sure ratio. 

In the study of the flow field within this flow regime, it is im- 
mediately recognized that an integral constraining relationship exists 
at the asymptotic state which points out the inadequacy of the indirect 
approach to the problem. Instead of asserting that the flow field is 
produced as a result of the upstream shock, it may be interpreted that 
the prevailing ambient and the asymptotic conditions forced the appear- 
ance of the shock with a specific geometry. This is, of course, typi- 
cal for all elliptic types of flow problems. 

The suggested method of solution to these problems relies on the 
selection of the asymptotic profile. It is not known whether there 
exist other approaches or methods which are equally effective in deal- 
ing with these problems. 
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Figure la Regular Reflection produced from Two-Dimensional 
Overexpanded Nozzle Flow 



Figure lb Mach Reflection produced from Two-Dimensional 
Overexpanded Nozzle Flow 



Figure lc Curved Shock produced from Two-Dimensional 
Overexpanded Nozzle Flow 
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Figure 4 Regions of Regular, Mach Reflection, and 

Curved Shock, p /p vs M [4] 
a oe 




Figure 6a Control Volume for Calculating the Asymptotic Velocity Profile 
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Figure 7 Shock Configuration at Corner of Nozzle 









Curved Shock Profile 


Asymptotic Velocity Profile 


Figure 9a Curved Shock and Asymptotic Velocity Profiles for M = 1.49 




ORIGINAL PAGE IS 
OF POOR QUALITY 







Figure 10b Variation of Centerline Pressure for = 3.0 and 5,0 
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Figure 13b Velocity u Distribution for M 
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Figure 15 Results of Velocity u Distribution 

along £ = 0.5 from Predictor and 

Final Calculations for M = 1.49 
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Table 1 C Values for Different Flow Conditions 


M = 
e 

1.49 

M = 
e 

2.0 

M = 
e 

3.0 

M = 

e 

4.0 

M = 
e 

5.0 

p /p 
1 a oe 

C 

P /P 
r a r oe 

C 

Pa^oe 

C 

P /p 
r oe 

C 

P /P 
r a r oe 

C 

0.50 

0.6067 

0.42 

1.4038 

0.22 

2.2994 

0.07 

1.0426 

0.044 

2.8406 

0.51 

0. 8410 

0.43 

1.5455 

0.23 

2.7216 

0.08 

1.5317 

0.045 

3.1077 

0.52 

1.0352 

0.44 

1.6923 

0.24 

3.3338 

0.09 

2.0884 

0.046 

3.4330 

0.53 

1.2121 

0.45 

1.8484 

0.25 

4.3267 

0.10 

2.9865 

0.047 

3.8388 

0.54 

1.3826 

0.46 

2.0191 

0.26 

6.2304 

0.11 

5.2914 

0.048 

4.3604 

0.55 

1.5540 

0.47 

2. 2111 

0.27 

11.4360 

0.12 

31.4710 

0.049 

5.0564 

0.56 

1.7331 

0.48 

2.4348 

0.28 

72.0330 

0.121 

56.700 

0.050 

6.0341 

0.57 

1.9276 

0.49 

2.7047 





0.051 

7.5120 

0.58 

2.1480 

0.50 

3.0427 





0.052 

10.0170 

0.59 

2.4097 

0.51 

3.4832 





0.053 

15.2500 

0.60 

2.7373 

0.52 

4.0826 





0.054 

32.9090 

0.61 

3.1709 

0.53 

4.9463 





0.0546 

54.6450 

0.62 

3.7792 

0.54 

6.2994 







0.63 

4.6953 

0.55 

8.7246 







0.64 

6. 2289 

0.56 

14.3560 







0.65 

9. 3159 

0.57 

40.3770 







0.66 

18.7750 

0.575 

75.8470 





- 


0.668 

77.1230 
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APPENDIX 

FLOW WITH STRONG CURVED SHOCK WITH AXIALLY SYMMETRIC CONFIGURATION 
The basic flow pattern of the axially symmetric flow within tills 
flow regimeis similar to that of two-dimensional flow. The curved 
shock occurs at the exit of the nozzle and becomes a normal shock at 
the axis of symmetry. Behind the shock, the flow is rotational and 
bounded by a free jet boundary. Eventually this free jet flow reaches 
an asymptotic state at far down stream position (x °°) where the 
pressure is constant. A schematic diagram depicting such a flow field 
is shown in Fig. 5 where x and y are interpreted as cylindrical coordi- 
nates . 


A, 1 BASIC RELATIONS GOVERNING THE FLOW 

The governing relations are listed sequentially in the following 
sections . 

A. 1.1 Shock Relations 

The governing shock equations are the same as those for 
two-dimensional flow. See Eqs. (2.1) through (2.5). 


A. 1.2 Free Jet Flow Region 

For axially symmetric inviscid flow of a compressible fluid 
the continuity and momentum principles can be given as: 

5puy . 9pvy _ 

3x 3y 


pu 
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pv 

3u 

1 3£ _ 

3x 

3y + 

.,2 3x 
K M 
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pu 

3v 
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3v 

1 j*£ = 
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.,2 3y 

if M J 
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(A.l) 

CA.2) 


0 


(A. 3) 
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The energy equation can be given as 

= o . (A. 4) 

Dt 

The flow region is bounded by (1) a constant pressure jet boundary, 

(2) the upstream curved shock, and (3) the asymptotic state prevailing 

at far downstream position (i.e,, x -*• °°) . 

The boundary condition for Eqs . (A.l), (A. 2), and [A- 3) are: 

p = p along the free jet boundary, 
a. 

u = u s , v = v s , p = P s » and P = P s along the curved shock, 
v = 0, |~ • = 0 at y = 0. 


A. 1.3 Constraining Condition at the Asymptotic State 

Consider the control volume shown in Fig. 19. The continuity 
principle would yield 


I 


y . 

aj 


/ 


p a u a y a dy a = 


CA.5) 


and the x-momentum principle also gives 


/ 


y . 
' aj 


2 j p a y e P e y e 2 >'e 

p a u a y a dy a + 7*72 = -~2 P e u e ~T 

2k M 2k M 
e e 


[A. 6) 


which can be rewritten as 


y aj 


r a J 2 . 1/, P a- ] \ 

I P * a * ir ' d, * , n l -T7) 


(A. 7) 


After introducing p u and y . as reference quantities for the 

aj aj aj 

integrands in Eqs. (A. 5) and (A. 7) and defining a » u /.u . and q = y/y 

a . aj ■ aj 
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these equations become 


p - u . y . 
aj aj aj 


— dn-j 
aj 


(A. 8) 


p - u . y . o 
aj a ^a 3 J p 

o 


$ I?) dp =-(! •- 


Pa “ 1 

.,2 

ic M 


(A. 9) 


Upon combining Eqs. (A. 8) and (A. 9), one obtains 

1 


1 - 


X 

I 


K M 


2 (P a - ^ 


P ’ 
aj 


4) n N> - 


dr) = 0 


(A. 10) 


K M 


o ' e 

A method similar to the two-dimensional flow method has been applied 
to calculate this flow case. Detailed calculation procedures will be 
discussed in the following sections. 


A. 2 ASYMPTOTIC STATE 

Since the solution of the problem relies on the establishment of 
the asymptotic flow profile, a detailed examination of the conditions 
imposed on such an asymptotic flow state is, therefore, necessary. 

Under the asymptotic flow condition, the slope of the velocity pro- 
file is precisely the vorticity of the flow. A study of the vorticity 
generated by the shock at the corner of the nozzle would also yield the 
slope of the profile of the jet boundary streamline at the asymptotic 
state. 


A. 2.1 Vorticity at the Jet Boundary Streamline 

Referring to Fig. 7, the continuity equation and the equation 
of motion in the streamline coordinate system are (in dimensional form): 
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1 - M 


e .90 _ sin 6 


Ss Bn 


3q 1 3p 

as p as 

2 30 9p 
pq 37 = ' 37 • 


y 


(A. 11) 
(A. 12) 
(A. 13) 


The vorticity can be expressed as 


_ 3q 90 _ p as 

W 3n ^ 9s pqR 9n 


(A. 14) 


For a j et boundary streamline at the corner of the nozzle, Eq. (A. 11) 
can be written as 

sin 9 


30 

an 


9£ = 
as 


(A, 15) 


At asymptotic state, the relation between the vorticity and the slope 
of velocity profile can be expressed as 


03 = 

a dn 


aj_ jchj) 


(A. 16) 


The vorticity at the shock can be written as 

■) ■ 


= _p_ 9S 
W s pqR 9n 


s p / as da ae' 

pqR V da d0 9n 


(A. 17) 


s s 

For axially symmetric flow, the vorticity along each streamline is 
not constant. It is proportional to the pressure moment y p relative 
to the axis of symmetry. Thus, along the constant pressure free jet 
boundary, the vorticity at the asymptotic state can be related to the 
vorticity at the shock by 


O) =03 

a 


^a _ p ( dS_ _dp _30\ ^a. _ 9q_ 
s y^ pqR \ da d0 9n j 9n 


(A. 18) 


The change of entropy with respect to the shock angle dS/da can be 
obtained from the following relations: 
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dS _ 4 Rk 

da (k - 1) (k + 1) 


M sin a cos a 
e 


2k %i 2 . 2 _ k - 1 

r* M sm a - - 7 ■ V 

k + 1 e k + 1 




cos a/(M 0 sin^ a) 


2 K - 1 

, , % . ,2 . 2 K + 1 

(k + 1) M sm a 


(A. 19) 


Taking the derivative of Eq. (2.2), one obtains 


da 

d6 


1/sin 2 6 


^ /— T* 

cos a 


K + 1 


M' 


2 ,,2 .2 . 
M sm a - 1 
e 


- 1 


(K + 1) sin 2 a 

? ? 2~ 
(M sin a - 1) 


(A. 20) 


thus, the. slope of the asymptotic velocity profile can be related to the 
shock angle, shock turning angle, and the asymptotic free jet height by 


M. 

3ri 


p sin 6 
r s 


4k 


= 1 " q ai + 13 


cos a 


M sin a cos a 
e 


.,2 • 3 „ 

M sm a 
e 


2k „,2 . 2 k - 1 

M sm a - 


k + 1 e 


K + 1 


2 K - 1 

. . .2 .2 K+l 

(k + 1) sm a 


- 2 - 
sm 6 


cos a 


K+l 


M‘ 


2 . 2 

2 M sin a - 1 
e 


- 1 


(k + 1) M 0 sin 2 a 
(M^ sin a - 1) 


2 

y . 

aj _ F 
2 

y„ 


F (M e , a, 6, y ) 


(A. 21) 
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It shou.ld be pointed out that at certain pressure ratios, da/d6 
in Eq. (A. 20) becomes undefined (i.e., da/d<5 = + 00 ) as shown in Fig. 20. 
This implies that the slope of the asymptotic velocity profile is also 
undefined. However, when the pressure ratios are different from that 
particular pressure ratio, present calculations are successful and the 
corresponding shock profiles can be computed. 


A. 2. 2 Velocity Profile at the Asymptotic State 

The velocity profile <j> in Eq. (A. 10) should satisfy the con- 
ditions of: 


(1) At q = 0: 


, , ac , n 

A = 6 = and ~r~ = 0 

T T c u . dp 

aj 


(2) At n = 1: 


* = 1 and It = F °V 6 ' y aj 3 • 


A similar velocity profile has been selected as 


<p - A + B erf [C (1 - q 2 )] + D erf [C (1 - q 4 ) ] . 


(A. 22) 


Since the boundary condition d<f)/dq 


^ is related to the asymptotic' 

free jet boundary height y iteration must be applied to determine the 

aj 

values of those coefficients and the free jet boundary height, y . . 

d I 


A. 2. 3 Corresponding Upstream Shock Profile 

After the asymptotic free jet height y . and its detailed pro- 

aj 

file have been established, the corresponding curved shock configuration 
can be determined through the following consideration. 



64 


If one carries out the integration, 


P a u a y dy 


along the asymptotic profile, continuity principle would imply that the 
fluid assuming the position y at the asymptotic section would have 

ci 

passed through the shock at height y s which is equal to 


y. 


1/2 


y = 2 
f s 


(/ u a y d * ) 


(A. 23) 


Since the velocity of the fluid at y has a known value of 4 , its stag- 

a a 

nation pressure p can be found from 


oa 


(k/k-1) 


i „2 .2 

l - c .4 
aj a 


(2.32) 


The ratio of the stagnation pressure across the shock at y s can be re- 
lated through the shock wave angle o by 


K + 1 .,2 . 2 

— - — M sm 0 
2 e 

, , K + 1 Xi 2 . 2 „ 
1 + — - — M sin a 
2 e 


(k/k-1) 


2k ,,2 . 2 . K - l\ 

M sm a - 

K+le k + 1 j 


P p P 

1 oa _ oa a 

(1/K-l) p P P 

1 oe 1 a oe 


(A. 24) 


The wave angle a, thus, can be found as a function of any position y 
along the shock. 

The calculated shock profiles and asymptotic velocity profiles for 
M e =1.49 at different pressure ratios are given in Fig. 20. 

It should be noted that at low pressure ratios, the maximum asymptotic 
velocity is not at free jet boundary. This is, of course, a result of the 
positive vorticity occurring at the corner of the nozzle under these specific 
pressure ratios. 
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